R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

# Lingting Shi Columbia Univeristy May 2019 How to run seruat
#Source https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html

#Clean working enviorment
rm(list=ls())
gc()

            used   (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
Ncells  10157230  542.5   15513551  828.6         NA  15513551  828.6
Vcells 278812045 2127.2  567139213 4327.0      65536 743937174 5675.8

#remotes::install_github('satijalab/seurat-wrappers')

#Install Package
#install.packages('Seurat')
#install.packages('dplyr')
#install.packages('umap')
#install.packages("reticulate")
#Load library
### Compare soft vs hard
#devtools::install_github('satijalab/seurat-data')
#install.packages("cowplot")
#install.packages("SCORPIUS")
#devtools::install_github('dviraran/SingleR')
#library(SingleR)
library(cowplot)
library(dplyr)
library(Seurat)
library(umap)
library(reticulate)
library(SCORPIUS)
library(ggplot2)
library(MASS) 
library(monocle3)
library(Seurat)
library(SeuratWrappers)
library(patchwork)

#devtools::install_github("cellgeni/sceasy")

# Crying, covert anndata to seurat
sceasy::convertFormat('/Users/lingtingshi/Documents/Azizi_lab/Decipher_monocle/gastric_data.h5', from="anndata", to="seurat",outFile='gastric_data.rds')

X -> counts

An object of class Seurat 
8705 features across 12632 samples within 1 assay 
Active assay: RNA (8705 features, 0 variable features)

#Reading in seurat object
gastric_data = readRDS('/Users/lingtingshi/Documents/Azizi_lab/Decipher_monocle/gastric_data.rds')

#Convert seurat object to cds
gastric_data@active.assay = 'RNA'
cds <- as.cell_data_set(gastric_data)

rowData(cds)$gene_name <- rownames(cds)
rowData(cds)$gene_short_name <- rowData(cds)$gene_name
rowData(cds)$id = rowData(cds)$gene_short_name
cds <- estimate_size_factors(cds)

cds<-preprocess_cds(cds, num_dim = 10)

cds <-reduce_dimension(cds, umap.min_dist = 1, umap.n_neighbors = 10L)

No preprocess_method specified, using preprocess_method = 'PCA'

cds <- cluster_cells(cds, resolution=1e-3)

plot_cells(cds, label_groups_by_cluster=FALSE,  color_cells_by = "class",group_label_size = 6, cell_size = 0.5)+
theme(legend.position = "right",
        legend.key.width = unit(0.2,"line"),
        legend.key.height = unit(0.2,"line"))

gastric_genes <- c("MUC6", "ENO1", "TFF3", "EPCAM", "MUC5AC",  "MUC13",
    "PGC", "ANXA2", "CEACAM5", "COL1A2", "TIMP1", "SPARC")
plot_cells(cds,
           genes=gastric_genes,
           show_trajectory_graph=FALSE)

plot_cells(cds, color_cells_by = "partition",graph_label_size = 50, cell_size = 0.5)

No trajectory to plot. Has learn_graph() been called yet?

cds <- learn_graph(cds)

  |                                                                                                                                                   
  |                                                                                                                                             |   0%
  |                                                                                                                                                   
  |=============================================================================================================================================| 100%

Warning: Argument `neimode' is deprecated; use `mode' instead

plot_cells(cds,
           color_cells_by = "class",
           label_groups_by_cluster=TRUE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           group_label_size= 4)+
theme(legend.position = "right",
        legend.key.width = unit(0.2,"line"),
        legend.key.height = unit(0.2,"line"))
ggsave("UMAP_traj_class.pdf", dpi = 600, width = 4.5, height = 3)

# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, class="Enteroendo"){
  cell_ids <- which(colData(cds)[, "class"] == class)
  
  closest_vertex <-
  cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-
  igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
  (which.max(table(closest_vertex[cell_ids,]))))]
  
  root_pr_nodes
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))

plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           graph_label_size=1.5)
ggsave("UMAP_by_pseudotime.pdf", dpi = 600, width = 4.5, height = 3)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.